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Abstract We numerically investigate the case of the pla¬ 
nar circular restricted three-body problem where the more 
massive primary is an oblate spheroid. A thorough numer¬ 
ical analysis takes place in the conhguration (x,y) and the 
(x, E) space in which we classify initial conditions of orbits 
into three categories: (i) bounded, (ii) escaping and (iii) col- 
lisional. Our results reveal that the oblateness coefficient has 
a huge impact on the character of orbits. Interpreting the col- 
lisional motion as leaking in the phase space we related our 
results to both chaotic scattering and the theory of leaking 
Hamiltonian systems. We successfully located the escape as 
well as the collisional basins and we managed to correlate 
them with the corresponding escape and collision times. We 
hope our contribution to be useful for a further understand¬ 
ing of the escape and collision properties of motion in this 
interesting version of the restricted three-body problem. 

Keywords Restricted three-body problem; Escape dynam¬ 
ics; Escape basins; Eractal basin boundaries 


1 Introduction 

One of most interesting topics in nonlinear dynamics is the 
issue of leaking or escaping orbits (e.g., Contopoulos (1990); 
Contopoulos & Kaufmann (1992); Contopoulos et al. (1993); 
Schneider et al. (2002); Altmann et al. (2013); Nagler et al. 
(2007)). Over the years many studies have been devoted on 
Hamiltonian systems with escapes (e.g., Barrio et al. (2009); 
Ernst & Peters (2014); Kandrup et al. (1999); Navarro & 
Henrard (2001); Zotos (2014a,b, 2015a)). However the topic 
of escaping orbits in Hamiltonian systems is by far less ex- 
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plored than the closely related problem of chaotic scatter¬ 
ing, where a test particle coming from inhnity approaches 
and then scatters off a complex potential (e.g., Bleher et al. 
(1989); Benet et al. (1996, 1998); Jung et al. (1999, 1995); 
Jung & Tel (1991); Seoane et al. (2006); Seoane & Sanjuan 
(2007)). 

Undoubtedly, one of the most important topics in dy¬ 
namical astronomy is the classical restricted three-body prob¬ 
lem (RTBP). The RTBP describes the motion of a test body 
with an inhnitesimal mass moving under the gravitational 
effects of a two primary bodies with hnite masses which 
move in circular orbits around their common center of mass 
(Szebehely, 1967). It is assumed that the third test body 
does not influence the motion of the two primary bodies. 
This problem has many applications in celestial mechanics, 
stellar systems, artihcial satellites, galactic dynamics, chaos 
theory and molecular physics and therefore is still a very 
active and stimulating held of research. 

Generally, the shapes of the two primary bodies in the 
classical version of the RTBP are assumed to be spherically 
symmetric. However, we found in nature that several celes¬ 
tial bodies, such as Saturn and Jupiter are sufficiently oblate 
(Beatty et al., 1999) or even triaxial. Eurthermore, irregu¬ 
lar shapes are also possible, especially in the case of mi¬ 
nor planets (e.g., Ceres) and meteoroids (Millis et al., 1987; 
Norton & Chitwood, 2008). The oblateness or triaxiality of 
a celestial body can produce perturbation deviations from 
the two-body motion. The study of oblateness coefficient 
includes the series of works of Beevi & Sharma (2012); 
Markellos et al. (1996, 2002); Kalantonis et al. (2005, 2006, 
2008); Kalvouridis & Gousidou-Koutita (2012); Perdiou et 
al. (2012); Sharma & Subba Rao (1979, 1986); Subba Rao 
& Sharma (1988, 1997); Sharma (1981, 1987, 1989, 1990); 
Singh &. Leke (2012, 2013) by considering the more mas¬ 
sive primary as an oblate spheroid with its equatorial plane 
co-incident with the plane of motion of the primaries. 
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All planets of the Solar System are in fact not spherically 
symmetric but oblate spheroidals. Consequently the exact 
shape of the planets and their oblateness must be taken into 
account when someone wants to model a space flight mis¬ 
sion. In the work of Oberti & Vienne (2003) it was proved 
that the theoretical data from the celestial systems Satum- 
Tethys-satellite and Saturn-Dione-satellite are much more 
accurate with respect with the corresponding observational 
data when the oblateness of Saturn is taken into considera¬ 
tion. Similarly, the oblateness of the Neptune was included 
in the work of Stuchi et al. (2008) regarding the dynamics 
of a spacecraft in the Neptune-Triton system. 

In the present paper we continue the work initiated in 
Zotos (2015b) following the same numerical techniques. Our 
aim is to numerically investigate how the oblateness coeffi¬ 
cient influences the nature of orbits in the restricted three- 
body problem. The paper is organized as follows: In Section 
2 we describe the properties of the considered dynamical 
model along with some necessary theoretical details. All the 
computational methods we used in order to determine the 
character of orbits are explained in Section 3. In the follow¬ 
ing Section, we conduct a thorough numerical investigation 
revealing the overall orbital structure (bounded regions and 
basins of escape/collision) of the system and how it is af¬ 
fected by the value of the oblateness coefficient with respect 
to the total orbital energy. Our paper ends with Section 5, 
where the main conclusions of this work are given. 

2 Details of the dynamical model 

It would be very illuminating to briefly recall the basic prop¬ 
erties and some aspects of the planar circular restricted three- 
body problem (Szebehely, 1967). The two primaries move 
on circular orbits with the same Kepler frequency around 
their common center of gravity, which is assumed to be fixed 
at the origin of the coordinates. The third body (test particle 
with mass much smaller than the masses of the primaries) 
moves in the same plane under the gravitational field of the 
two primaries. The non-dimensional masses of the two pri¬ 
maries are 1 - ju and ju, where n — m 2 l{m\ + m 2 ), with 
m\ > m 2 is the mass ratio. We consider an intermediate 
value of the mass ratio, that is ju = 1/10. This value remains 
constant throughout the paper. 

We choose as a reference frame a rotating coordinate 
system where the origin is at (0,0), while the centers C\ 
and C 2 of the two primary bodies are located at i-fi, 0 ) and 
(1 - /i, 0), respectively. The total time-independent effective 
potential according to Shatma & Subba Rao (1976) is 

s (l-A*) (l-/t)Ai n^(2^ 2\ 

V{x,y)^ - -3 - —[x^+y], (1) 

r 2 ri 2rf 2 '' 
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are the distances to the respective primaries and the angular 
velocity (n), while Ai is the oblateness coefficient which is 
defined as 
^ (R£)2 - (RP)^ 

^ 5R2 ’ 

where RE and RP are the equatorial and polar radius, re¬ 
spectively of the oblate primary, while R is the distance be¬ 
tween the centers of the two primaries. We consider values 
of the oblateness coefficient is in the interval [ 0 , 0 . 1 ] (see 
e.g., (Kalantonis et al., 2006; Perdios & Kalantonis, 2006)), 
while we study the effect of oblateness up to the linear coef¬ 
ficient J 2 only (Abouelmagd, 2012). 

The human species stands on the edge of a new fron¬ 
tier, the transition from a planet-bound to a space-faring civ¬ 
ilization. The expansion into the Solar System, in terms of 
dynamics of artificial satellites, requires the formulation of 
new models that include the effects of some of the perturb¬ 
ing forces on the satellite. Therefore, the problem is not only 
of mathematical interest but has astrophysical applications. 
The most striking example of perturbations arising from the 
oblateness in the Solar System is the orbit of the fifth satel¬ 
lite of Jupiter, Amalthea. This planet is so oblate and the 
satellite’s orbit is so small that its line of apsides advances 
about 900° in a year (e.g., Moulton (1914)). It should be 
pointed out that in general terms the values of the oblate¬ 
ness in the Solar System are relatively low (see e.g., Milone 
& Wilson (2014) for updated values of the oblateness co¬ 
efficient). There are however close-in extrasolar planetary 
systems with hot Jupiters in which the oblateness coefficient 
varies in the interval [0.001, 0.01]. In this work we consider 
even higher values of the oblateness coefficient (Ai < 0 . 1 ) 
because in the near future it is possible new highly oblate 
extrasolar planets to be discovered. 

The scaled equations of motion describing the motion of 
the test body in the corotating frame read (Sharma & Subba 
Rao, 1976) 


X - 2ny — 


dV(x,y) 

dx 


y — —2nx— 


dV(x,y) 

dy 


(4) 


The dynamical system (4) admits the well know Jacobi in¬ 
tegral 

J(x, y, x,y)^^ {x^ + f) + V(x, y) = E, (5) 

where x and y are the velocities, while E is the numerical 
value of the orbital energy which is conserved and defines a 
three-dimensional invariant manifold in the total four-dimensional 
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phase space. Thus, an orbit with a given value of it’s en¬ 
ergy integral is restricted in its motion to regions in which 
E < V{x,y), while all other regions are forbidden to the test 
body. It is widely believed that J is the only independent in¬ 
tegral of motion for the RTBP system (Poincare, 1993). The 
value of the energy E is related with the Jacobi constant by 
C = -2E. 

The dynamical system has five equilibria known as La- 
grangian points (Szebehely, 1967) at which 
dV{x,y)_dV{x,y)_ 

dx dy ' ^ ^ 

Three of them, Li, Li, and L 3 , are collinear points located 
in the jr-axis. The central stationary point Li is a local min¬ 
imum of the potential V(x, y). The stationary points L 2 and 
L 3 are saddle points. Let L 2 located at x > 0, while L 3 be 
at X < 0. The points L 4 and Lg on the other hand, are lo¬ 
cal maxima of the gravitational potential, enclosed by the 
banana-shaped isolines. The Lagrangian points are very im¬ 
portant especially for astronautical applications. This can be 
seen in the Sun-Jupiter system where several thousand as¬ 
teroids (collectively referred to as Trojan asteroids), are in 
orbits of equilibrium points. In practice these Lagrangian 
points have proven to be very useful indeed since a space¬ 
craft can be made to execute a small orbit about one of these 
equilibrium points with a very small expenditure of energy 
(Singh & Leke, 2014). 

The projection of the four-dimensional phase space onto 
the configuration (or position) space (x,y) is called the Hill’s 
regions and is divided into three domains: (i) the interior 
region for x(L 3 ) < x < x(L 2 ); (ii) the exterior region for 
X < x(L 3 ) and x > x(L 2 ); (iii) the forbidden regions. The 
boundaries of these Hill’s regions are called Zero Velocity 
Curves (ZVCs) because they are the locus in the configura¬ 
tion (x,y) space where the kinetic energy vanishes. 

The values of the Jacobi integral at the five Lagrangian 
points Li, / = 1, ...,5 are critical energy levels and they are 
denoted as Ei (Note that £4 = £ 5 ). The structure of the 
equipotential surfaces strongly depends on the value of the 
energy. In particular, there are five distinct cases 

- £ < £ 1 : All necks are closed, so there are only bounded 
and collisional basins. 

-El < £ < £ 2 : Only the neck around £1 is open thus 
allowing orbits to move around both primaries. 

- £2 < £ < £3: The neck around £2 is open, so orbits can 
enter the exterior region and escape form the system. 

- £3 < £ < £ 4 : The necks around both £2 and £3 are 
open, so orbits can escape through two different escape 
channels. 

- £ > £ 4 : The banana-shaped forbidden regions disappear 
and therefore, motion over the entire configuration (x, y) 
space is possible. 

In Fig. l(a-d) we present for Ai - 0.01 the structure of the 
first four possible Hill’s region configurations. We observe 


in Fig. Id the two openings (exit channels) at the Lagrangian 
points £2 and £3 through which the test body can enter the 
exterior region and then leak out. In fact, we may say that 
these two exits act as hoses connecting the interior region 
of the system where x(£ 3 ) < x < x(£ 2 ) with the “out¬ 
side world” of the exterior region. The position of the La¬ 
grangian points as well as the critical values of the energy 
are functions of the oblateness coefficient Ai (e.g., Singh 
& Leke (2012)). In Table 1 we provide the location of the 
Lagrangian points and the critical values of the total orbital 
energy when Ai = {0,0.001,0.01,0.1). 

3 Computational methods and criteria 

The motion of the test third body is restricted to a three- 
dimensional surface £ = const, due to the existence of the 
Jacobi integral. With polar coordinates (r, (p) in the center of 
the mass system of the corotating frame the condition r - 0 
defines a two-dimensional surface of section, with two dis¬ 
joint parts (^ < 0 and 0 > 0. Each of these two parts has a 
unique projection onto the configuration (x,y) space. In or¬ 
der to explore the orbital structure of the system we need to 
define samples of initial conditions of orbits whose proper¬ 
ties will be identified. Fore this purpose, we define for sev¬ 
eral values of the total orbital energy £, as well as for the 
oblateness coefficient Ai dense uniform grids of 1024x 1024 
initial conditions regularly distributed on the configuration 
(x, y) space inside the area allowed by the value of the to¬ 
tal orbital energy. Following a typical approach, the orbits 
are launched with initial conditions inside a certain region, 
called scattering region, which in our case is a square grid 
with -2 < x,y < 2 . 

In the RTBP system the configuration space extends to 
infinity thus making the identification of the type of mo¬ 
tion of the test body for specific initial conditions a rather 
demanding task. There are three possible types of motion 
for the test body: (i) bounded motion around one of the 
primaries, or even around both; (ii) escape to infinity; (iii) 
collision into one of the two primaries. Now we need to 
define appropriate numerical criteria for distinguishing be¬ 
tween these three types of motion. The motion is considered 
as bounded if the test body stays confined for integration 
time fmax inside the system’s disk with radius Rd and cen¬ 
ter coinciding with the center of mass origin at (0,0). Ob¬ 
viously, the higher the values of fmax and Rd the more plau¬ 
sible becomes the definition of bounded motion and in the 
limit fmax ^ °° the definition is the precise description of 
bounded motion in a finite disk of radius Rd- Consequently, 
the higher these two values, the longer the numerical inte¬ 
gration of initial conditions of orbits lasts. In our calcula¬ 
tions we choose fmax = 10"^ and Rd - 10 (see Fig. 2) as in 
Nagler (2004, 2005) and Zotos (2015b). We decided to in¬ 
clude a relatively high disk radius {Rd - 10) in order to be 
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Fig. 1 Four possible Hill’s regions configurations for the RTBP system when Ai = 0.01. The white domains correspond to the Hill’s regions, 
gray shaded domains indicate the forbidden regions, while the thick black lines depict the Zero Velocity Curves (ZVCs). The red dots pinpoint the 
position of the Lagrangian points, while the positions of the centers of the two primary bodies are indicated by blue dots, (a-upper left): E = -1.92; 
(b-upper right): E = -1.78; (c-lower left): E = -1.73; (d-lower right): E = -1.55. 


sure that the orbits will certainly escape from the system and 
not return back to the interior region. Furthermore, it should 
be emphasized that for low values of f^ax the fractal bound¬ 
aries of stability islands corresponding to bounded motion 
become more smooth. Moreover, an orbit is identified as es¬ 
caping and the numerical integration stops if the test body 
body intersects the system’s disk with velocity pointing out¬ 
wards at a time fesc < Wx- Finally, a collision with one of 
the primaries occurs if the test body, assuming it is a point 


mass, crosses the disk with radius around the primary. 
For the larger oblate primary we choose - 10 Gener¬ 
ally, it is assumed that the radius of a celestial body (e.g., a 
planet) is directly proportional to the cubic root of its mass. 
Therefore, for the sake of simplicity of the numerical calcu¬ 
lations we adopt the simple relation between the radii of the 
primaries 


Rm, = Rnn X (2^)h3 , 


(7) 
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Table 1 The position of the Lagrangian points and the critical values of the total orbital energy for four values of the oblateness coefficient. 


Al 

Al 

^2 

Li 

u 

0.000 

(0.60903511,0) 

(1.25969983,0) 

(-1.04160890, 0) 

(0.40000000, 0.86602540) 

0.001 

(0.60934671, 0) 

(1.25945159, 0) 

(-1.04165743,0) 

(0.40049937, 0.86573689) 

0.010 

(0.61207238, 0) 

(1.25724426, 0) 

(-1.04208268,0) 

(0.40493832, 0.86315542) 

0.100 

(0.63363978, 0) 

(1.23751047,0) 

(-1.04544100, 0) 

(0.44448280, 0.83877200) 

Al 

El 

E2 

Ei 

£4 

0.000 

-1.79847661 

-1.73334221 

-1.54978907 

-1.45500000 

0.001 

-1.80001655 

-1.73471117 

-1.55114179 

-1.45613246 

0.010 

-1.81381551 

-1.74701526 

-1.56331597 

-1.46632128 

0.100 

-1.94701666 

-1.86855958 

-1.68503212 

-1.56790343 


Outside the system’s disk: escape orbit 



Fig. 2 Schematic picture of the three different types of motion. The 
motion is considered to be bounded if the test body stays confined for 
integration time inside the system’s disk with radius Rj = 10, 
while the motion is unbounded and the numerical integration stops 
when the test body crosses the system’s disk with velocity pointing 
outwards. Collision or (crash) with one of the primaries occurs when 
the test body crosses the disk of radii R,„i and R „2 of the primaries. 


which was introduced in Nagler (2005). In Nagler (2004, 
2005) it was shown that the radii of the primaries influence 
the area of collision and escape basins. 

As it was stated earlier, in our computations, we set 10"^ 
time units as a maximum time of numerical integration. The 
vast majority of escaping orbits (regular and chaotic) how¬ 
ever, need considerable less time to escape from the system 
(obviously, the numerical integration is effectively ended when 
an orbit moves outside the system’s disk and escapes). Nev¬ 
ertheless, we decided to use such a vast integration time just 
to be sure that all orbits have enough time in order to es¬ 
cape. Remember, that there are the so called “sticky orbits” 
which behave as regular ones during long periods of time. 
Here we should clarify, that orbits which do not escape after 


a numerical integration of 10“* time units are considered as 
non-escaping or trapped. 

The equations of motion (4) for the initial conditions 
of all orbits are forwarded integrated using a double pre¬ 
cision Bulirsch-Stoer FORTRAN 77 algorithm (e.g.. Press et 
al. (1992)) with a small time step of order of 10“^, which 
is sufficient enough for the desired accuracy of our compu¬ 
tations. Here we should emphasize, that our previous nu¬ 
merical experience suggests that the Bulirsch-Stoer integra¬ 
tor is both faster and more accurate than a double precision 
Runge-Kutta-Fehlberg algorithm of order 7 with Cash-Karp 
coefficients. Throughout all our computations, the Jacobian 
energy integral (Eq. (5)) was conserved better than one part 
in 10“^', although for most orbits it was better than one part 
in 10 For collisional orbits where the test body moves 
inside a region of radius 10“^ around one of the primaries 
the Lemaitre’s global regularization method is applied. 


4 Numerical results & Orbit classification 

The main numerical task is to classify initial conditions of 
orbits in the <^ < 0 part' of the surface of section r - 0 into 
three categories: (i) bounded orbits; (ii) escaping orbits and 
(iii) collisional orbits. Moreover, two additional properties 
of the orbits will be examined: (i) the time-scale of collision 
and (ii) the time-scale of the escapes (we shall also use the 
terms escape period or escape rates). In this work we shall 
explore these dynamical quantities for various values of the 
total orbital energy, as well as for the oblateness coefficient 
Aj. In particular, we choose three energy levels which cor¬ 
respond to the last three Hill’s regions configurations. The 
first two Hill’s regions configurations contain only bounded 
and collisional orbits, so the orbital content is not so inter¬ 
esting. In the following color-coded grids (or orbit type di¬ 
agrams - OTDs) each pixel is assigned a color according to 
the orbit type. Thus the initial conditions of orbits are classi¬ 
fied into bounded orbits, unbounded or escaping orbits and 
collisional orbits. In this special type of Poincare surface of 

* We choose the 0 < 0 instead of the 0 > 0 part simply because in 
Zotos (2015b) we seen that it contains more interesting orbital content. 
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section the phase space emerges as a close and compact mix 
of escape basins, collisional basins and stability islands. Our 
numerical calculations indicate that apart from the escaping 
and collisional orbits there is also a considerable amount of 
non-escaping orbits. In general terms, the majority of non¬ 
escaping regions corresponds to initial conditions of regular 
orbits, where an adelphic integral of motion is present, re¬ 
stricting their accessible phase space and therefore hinders 
their escape. 

4.1 Results for Region III: E 2 < E < 

Our numerical exploration begins in the energy region III 
(£2 < E < £ 3 ) and we choose the energy level £ = -1.70. 
In Fig. 3 the OTD decompositions of the 0 < 0 part of the 
surface of section r - 0 reveal the orbital structure of the 
conhguration (x,y) space for four values of the oblateness 
coefficient. The black solid lines denote the ZVC, while the 
inaccessible forbidden regions are marked in gray. The color 
of a point represents the orbit type of a test body which has 
been launched with pericenter position at (x,y). When the 
oblateness coefficient is zero (Ai) = 0 we observe in Fig. 3a 
the followings: (i) around the centers of the two primaries 
there are stability islands, while in the exterior region an¬ 
other ring-shaped stability island is present. The stability 
islands located in the interior region correspond to regular 
orbits around one of the primary bodies, while the stability 
island in the exterior region is formed by initial conditions 
of regular orbits that circulate around both primaries; (ii) 
The region in the boundaries of the stability islands of the 
interior region is mainly occupied by orbits which lead to 
collision; (iii) Near the center of primary 1 we can identify 
a small hole which contains a mixture of collisional and es¬ 
caping orbits; (iv) outside the forbidden regions there is a 
well-formed circular basin of escaping orbits. It is evident 
in Fig. 3b that even a very low values of the oblateness coef- 
hcient (Ai = 0 . 001 ) is able to influence the orbital content. 
In particular, there are two main differences with respect to 
the case of zero oblateness shown in Fig. 3a. First the area 
of the stability region around the center of the oblate pri¬ 
mary is smaller, while the central hole contains now only 
initial conditions of collisional orbits. Moreover, the colli¬ 
sional basin around primary 1 is larger in size. The orbital 
structure changes even further when Ai = 0.01 as we can 
see in Fig. 3c. Almost all the area around the oblate primary 
1 is occupied by initial conditions of collisional orbits, while 
the corresponding stability island splits into two small parts. 
It seems that so far the region around smaller primary 2 as 
well as the exterior region remain almost unperturbed by the 
shift on the value of the oblateness coefficient. Things how¬ 
ever change drastically when the larger primary is highly 
oblate with Ai = 0.1. It is seen in Fig. 3d that the entire in¬ 
terior region around primary 1 is dominated by a large colli¬ 


sional basin. The exterior region on the other hand, is highly 
fractal, while the stability island of regular orbits circulating 
around both primaries is now absent. At least at this energy 
level our computations suggest that the increase in the value 
of the oblateness coefficient does not affect the stability is¬ 
land around primary 2 , however at the highest studied value 
the collisional basin around this stability island is weakened. 
It should be noted that when Ai = 0.1 several spiral colli¬ 
sional basins are inside the fractal exterior region. 

The following Fig. 4 shows how the escape and colli¬ 
sional times of orbits are distributed on the configuration 
{x,y) space for the four values of the oblateness coefficient 
discussed in Fig. 3. Light reddish colors correspond to fast 
escaping/collional orbits, dark blue/purple colors indicate 
large escape/collional times, while white color denote sta¬ 
bility islands of regular motion. Note that the scale on the 
color bar is logarithmic. Inspecting the spatial distribution 
of various different ranges of escape time, we are able to 
associate medium escape time with the stable manifold of 
a non-attracting chaotic invariant set, which is spread out 
throughout this region of the chaotic sea, while the largest 
escape time values on the other hand, are linked with sticky 
motion around the stability islands of the two primary bod¬ 
ies. As for the collisional time we see that orbits with ini¬ 
tial conditions very close to the vicinity of the center of the 
oblate primary 1 collide with it almost immediately, within 
the first time step of the numerical integration, while this 
phenomenon is not observed for the case of primary body 2 
which is not oblate (see also Zotos (2015b)). Looking more 
carefully Fig. 4d we observe that when Ai = 0.1 the area 
of the stability region around primary 2 (the only stability 
region that survives) is smaller with respect to the tree pre¬ 
vious cases shown in Figs. 4(a-c). Thus we may say that 
high enough values of the oblateness coefficient influence 
also the stability region corresponding to smaller spherically 
symmetric primary. 

4.2 Results for Region IV: £3 < £ < £4 

We continue our investigation in the energy region IV (£3 < 
£ < £ 4 ). It is evident from the critical energy values given 
in Table 1 that there is not a single energy level which we 
can use for all four cases regarding the value of the oblate¬ 
ness coefficient. Therefore for the first three cases, that is 
Ai = {0,0.001,0.01), we choose the value £1 = -1.54, 
while for the last case (Ai = 0 . 1 ) we use the energy level 
£1 = -1.67. The orbital structure of the conhguration (x,y) 
space is unveiled in Fig. 5 through the OTD decompositions 
of the (^ < 0 part of the surface of section r - 0. We observe 
in Fig. 5a where Ai = 0 that in this case the stability island 
in the exterior region is no longer present. The vast major¬ 
ity of the exterior region is occupied by initial conditions of 
escaping orbits, while inside the extended escape basin we 
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X X 

Fig. 3 The orbital structure of the ip < Q part of the surface of section r = 0 when £ = -1.7. (a-upper left): Ai = 0; (b-upper right): Ai = 0.001; 
(c-lower left): Ai = 0.01; (d-lower right): Ai =0.1. The color code is the following: bounded orbits (green), collisional orbits to oblate primary 1 
(blue), collisional orbits to primary 2 (red), and escaping orbits (cyan). 


identify delocalized initial conditions of orbits that collide to 
the oblate primary. On the other hand, the initial conditions 
of orbits which collide to primary 2 form well-defined spiral 
basins. The interior region is almost the same with that dis¬ 
cussed earlier in Fig. 3a, however there is a minor difference; 
the collisional basin around the stability island of primary 2 
is absent now. With the introduction of the oblateness coef¬ 


ficient it is seen in Figs. 5(b-d) that the interior region dis¬ 
plays exactly the same behaviour as in the previous case. In 
the exterior region two phenomena take place as the value 
of the oblateness coefficient increases: (i) the amount of the 
escaping orbits decreases while at the same time the portion 
of the collisional orbits to oblate primary 1 increases and 
(ii) the fractality of the exterior region gradually decreases 
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(a) (b) 



Fig. 4 Distribution of the escape and collisional time of the orbits on the <^ < 0 part of the surface of section r = 0 when E = -1.70 for the values 
of the oblateness coefficient of Fig. 3. The darker the color, the larger the escape/collisional time. Initial conditions of bounded regular orbits are 
shown in white. 


with increasing Ai and several escape and collisional basins 
emerge. 

The distribution of the escape and collisional times of 
orbits on the conhguration space is shown in Fig. 6 . One 
may observe that the results are very similar to those pre¬ 
sented earlier in Fig. 4, where we found that orbits with ini¬ 
tial conditions inside the escape and collisional basins have 
the smallest escape/collision rates, while on the other hand, 
the longest escape/collisional rates correspond to orbits with 
initial conditions in the fractal regions of the OTDs. Our cal¬ 
culations reveal, and this can be seen better in Figs. 6 (a-d), 
that in this energy region the oblateness coefficient does not 
affect the size of the stability island of motion around pri¬ 
mary 2 . 


4.3 Results for Region V; £ > £4 

The last case under consideration involves the results in the 
energy region V where the entire configuration space is avail¬ 
able for motion since the forbidden regions completely dis¬ 
appear. Once more, all the different aspects of the numeri¬ 
cal approach remain exactly the same as in the two previ¬ 
ously studied cases. Fig. 7 presents the orbital structure of 
the configuration space through the OTD decompositions of 
the Ip < 0 part of the surface of section r - 0. When both 
primaries are spherically symmetric, that is when Ai = 0 , 
we see in Fig. 7a that most of the configuration {x,y) plane 
is highly fractal, while a basin of collisional orbits is located 
around the stability island of primary 1. Around primary 2 
there is another stability island but now the corresponding 
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Fig. 5 The orbital structure of the 0 < 0 part of the surface of section f = 0 when (a-c): £ = -1.54 and (d-lower right): Ei = -1.67. (a-upper left): 
A[ = 0; (b-upper right): Ai = 0.001; (c-lower left): Ai = 0.01; (d-lower right): A[ =0.1. The color code is the same as in Fig. 3. 


collisional basin identified in energy region III is not present 
now. In addition, the central hole inside the stability island 
of primary 1 still survives. At the lowest value of the oblate¬ 
ness coefficient (Ai = 0.001) the observed changes in Fig. 
7b are the same with those reported earlier in Fig. 5b. The 
same applies in the case for Ai = 0.01 where the results 
are shown in Fig. 7c. Finally in Fig. 7d where we have the 


scenario with the highest possible value of the oblateness 
coefficient (A; = 0.1) it is seen that the vast majority of the 
(x,y) plane is occupied by well-formed escape basins and 
collisional basins, while the fractal areas are confined to the 
boundaries between the several basins. It should be pointed 
out that in this case the area of the stability region around 
primary 2 is considerably reduced. Thus one may reasonably 
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X X 

Fig. 6 Distribution of the escape and collisional time of the orbits on the 0 < 0 part of the surface of section r = 0 when (a-c): E = -1.54 and 
(d-lower right): £i = -1.67 for the values of the oblateness coefficient of Fig. 5. 


conclude that in the energy region V and for high enough 
values of Ai the stability region around spherically primary 
2 is highly affected by the oblateness coefficient. 

In Fig. 8 we depict the distribution of the escape and col¬ 
lisional times of orbits on the configuration space. One can 
see similar outcomes with that presented in the two previ¬ 
ous subsections. At this point, we would like to emphasize 
that the basins of escape can be easily distinguished in Fig. 
8 , being the regions with intermediate colors indicating fast 
escaping orbits. Indeed, our numerical computations suggest 
that orbits with initial conditions inside these basins need no 
more than 10 time units in order to escape from the system. 
Furthermore, the collisional basins are shown with reddish 
colors where the corresponding collisional time is less than 
one time unit of numerical integration. 


The OTDs shown in Figs. 3, 5 and 7 have both frac¬ 
tal and non-fractal (smooth) boundary regions which sepa¬ 
rate the escape basins from the collisional basins. Such frac¬ 
tal basin boundaries is a common phenomenon in leaking 
Hamiltonian systems (e.g., Bleher et al. (1988); de Moura & 
Letelier (1999); de Moura & Grebogi (2002); Schneider et 
al. (2002, 2003); Tuval et al. (2004)). In the RTBP system 
the leakages are defined by both escape and collision condi¬ 
tions thus resulting in three exit modes. However, due to the 
high complexity of the basin boundaries, it is very difficult, 
or even impossible, to predict in these regions whether the 
test body (e.g., a satellite, asteroid, planet etc) collides with 
one of the primary bodies or escapes from the dynamical 
system. 
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X X 

Fig. 7 The orbital structure of the 0 < 0 part of the surface of section r = 0 when £i =-1.45. (a-upper left): Ai = 0; (b-upper right): Ai =0.001; 
(c-lower left): Ai = 0.01; (d-lower right): Ai =0.1. The color code is the same as in Fig. 3. 


4.4 An overview analysis 

The color-coded OTDs in the configuration (x, y) space pro¬ 
vide sufficient information on the phase space mixing how¬ 
ever, for only a fixed value of the energy integral and also 
for orbits that traverse the surface of section retrogradely. 
Henon (Henon, 1969), introduced a new type of plane which 
can provide information not only about stability and chaotic 


regions but also about areas of bounded and unbounded mo¬ 
tion using the section y = i = 0, y>0 (see also Barrio 
et al. (2008)). In other words, all the initial conditions of 
the orbits of the test particles are launched from the ;ic-axis 
with X = xq, parallel to the y-axis (y = 0). Consequently, 
in contrast to the previously discussed types of planes, only 
orbits with pericenters on the .r-axis are included and there- 
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X X 

Fig. 8 Distribution of the escape and collisional time of the orbits on the (p <0 part of the surface of section r = 0 when E\ = -1.45 for the values 
of the oblateness coefficient of Eig. 7. 


fore, the value of the energy E can now be used as an ordi¬ 
nate. In this way, we can monitor how the energy influences 
the overall orbital structure of our dynamical system using 
a continuous spectrum of energy values rather than few dis¬ 
crete energy levels. In Fig. 9 we present the orbital structure 
of the (x, E) plane for four values of the oblateness coeffi¬ 
cient when E € [-3,1], while in Fig. 10 the distribution of 
the corresponding escape and collision times of the orbits is 
depicted. 

We can observe the presence of several types of regu¬ 
lar orbits around the two primary bodies. Being more pre¬ 
cise, on both sides of the primaries we identify stability is¬ 
lands corresponding to both direct (counterclockwise) and 
retrograde (clockwise) quasi-periodic orbits. It is seen that 
a large portion of the exterior region, that is for x < x(Lj) 
and X > x(L 2 ), a large portion of the (x, E) plane is covered 


by initial conditions of escaping orbits however, at the left- 
hand side of the same plane two stability islands of regular 
orbits that circulate around both primaries are observed. Ad¬ 
ditional numerical calculations reveal that for much lower 
values of X (x < 5) these two stability islands are joined 
and form a crescent-like shape. We also see that collisional 
basins to oblate primary 1 leak outside the interior region, 
mainly outside L 2 , and create complicated spiral shapes in 
the exterior region. On the other hand, the thin bands rep¬ 
resent initial conditions of orbits that collide with primary 
body 2 are much more confined. It should be pointed out 
that in the blow-ups of the diagram several additional very 
small islands of stability have been identified^. 


^ An infinite number of regions of (stable) quasi-periodic (or small 
scale chaotic) motion is expected from classical chaos theory. 
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As the value of the oblateness coefficient increases the 
structure of the {x, E) planes exhibits the following changes: 
(i) The collisional basins to oblate primary 1 significantly 
grow and for Ai = 0.1 they dominate; (ii) The area of the 
stability islands around primary 1 reduces and at the highest 
studied value of Ai they completely disappear. According 
to Broucke’s classification (Broucke, 1968) the periodic or¬ 


bits around the primaries belong to the families C (at the 
left side of the primary) and Hi (at the right side of the 
primary), while Simo & Stuchi (2000) proved for the pla¬ 
nar Hill’s problem that the stability regions of the C family 
are more stable than those of the Hi family. Our computa¬ 
tions show that both families disappear for large values of 
the oblateness coefficient; (iii) The stability islands around 
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Fig. 10 Distribution of the escape and collisional time of the orbits on the (x, E) plane for the values of the oblateness coefficient of Fig. 9. 


primary body 2 are almost unperturbed, at least in the inter¬ 
val E € [-3, -1.5], by the shifting on the value of the oblate¬ 
ness. The phenomenon that stability islands can appear and 
disappear as a dynamical parameter is changed has also been 
reported in earlier paper (e.g.. Barrio et al. (2006); de Assis 
&. Terra (2014)); (iv) Inspecting Figs. lO(a-d) we see that 
orbits with initial conditions very close to the center of the 
oblate primary 1 collide almost immediately with it, while 
their portion (or in other words the thickness of the verti¬ 
cal basin) increases for larger values of the oblateness co¬ 
efficient Aj; (v) Another interesting phenomenon is the fact 
that as the primary body 1 becomes more and more oblate 
the fractality of the (x, E) plane reduces and the boundaries 
between escaping and collisional basins appear to become 
smoother. It should be emphasized that the fractality of the 
structures was not measured by computing the correspond¬ 


ing fractal dimension. When we state that an area is fractal 
we mean that it has a fractal-like geometry. 

It would be very informative to monitor the evolution of 
the percentages of all the different types of orbits as a func¬ 
tion of the total orbital energy E for the (x, E) planes shown 
in Figs. 9(a-d). Our results are presented in Figs, ll(a-d). 
We see that in all four cases the rate of the escaping orbits 
starts at about 80% for E - -3 and gradually reduces until 
about 20% for E - -1.8. For E > -1.8 it suddenly in¬ 
creases but this increase is controlled by the oblateness co¬ 
efficient. In particular, the higher the value of A i the lower 
the increase of the escaping percentage in the energy inter¬ 
val [-1.8,1]. In the same interval however, according to Fig. 
10a escaping orbits dominate the (x, E) plane when Ai = 0. 
The evolution of the percentage of regular orbits displays 
similar patterns in the energy interval [-3,-1.8] however, 
the corresponding rates are reduced with increasing Ai. For 






































Influence of the oblateness coefficient in the restricted three-body problem 


15 





Fig. 11 Evolution of the percentages of escaping, regular and collisional orbits on the (x, £)-plane as a function of the total orbital energy E. 
(a-upper left): Ai = 0; (b-upper right): Aj = 0.001; (c-lower left): Ai = 0.01; (d-lower right): Ai =0.1. 


larger values of the energy the percentage of bounded or¬ 
bits fluctuate and in general terms we may say that again 
the oblateness coefficient has an influence. The percentage 
of collisional orbits to oblate primary 1 has a monotone be¬ 
haviour in the energy interval [-3,-1.8], while for £ > -1.8 
it increases. This increase becomes stronger and stronger as 
the value of the oblateness coefficient increases. For large 


enough values of the energy the rate of collisional orbits to 
primary 1 decreases in all four cases, while at the same time 
escaping orbits is the most populated type of orbits. The per¬ 
centage of collisional orbits to primary 2 on the other hand, 
is extremely low, while it peaks only in the energy interval 
[-1.8, -1.2]. Once more, the peaks in this energy range are 
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(a) 




E 


Fig. 12 (a-left): The average escape time of orbits < t^sc > and (b-right): the average collision time of orbits < tcoi > as a function of the total 
orbital energy E. 


influenced by the value of Ai; the larger the oblateness co- 
efflcient the lower the peaks. 

Another interesting aspect would be to reveal how the 
oblateness coefficient influences the escape as well as the 
collision time of orbits. The evolution of the average value 
of the escape time < fgsc > of orbits as a function of the total 
orbital energy is given in Fig. 12a. It is evident, especially in 
the energy interval [-1.8,1], that as the value of the oblate¬ 
ness coefficient increases the escape rates of orbits are re¬ 
duced. In the same vein in Fig. 12b we present the evolution 
of the average collision time < fcoi > of orbits that collide to 
oblate primary 1. Once more we observe that the collision 
time of orbits decreases with increasing value of Ai, espe¬ 
cially in the energy interval [-1.6,1]. Additional numerical 
calculations, not shown here, indicate that the oblateness co¬ 
efficient also influences the collision rates of orbits which 
collide to spherical primary 2. In particular, there are energy 
ranges in which the collision rates to primary 2 are reduced 
with increasing A 1 , however the diagrams are not so clear as 
in Fig. 12b, so we decided not to include them. 

Before closing this section we would like to add that the 
particular value of the mass ratio fj. does not really change 
the qualitative nature of the numerical outcomes presented 
in this section. Indeed after conducting some additional cal¬ 
culations with larger and lower values of fj. we concluded 
to the same results. The parameters which mostly influence 
the orbital dynamics are the total energy and of course the 
oblateness coefficient. 


5 Conclusions 

The aim of this work was to reveal how the oblateness coeffi¬ 
cient influences the character of orbits in the classical planar 
circular restricted three-body problem. After conducting an 
extensive and thorough numerical investigation we managed 
to distinguish between bounded, escaping and collisional or¬ 
bits and we also located the basins of escape and collision, 
finding also correlations with the corresponding escape and 
collision times. Our numerical results strongly suggest that 
the oblateness coefficient plays a very important role in the 
nature of the test’s body motion under the gravitational field 
of two primaries. To our knowledge, this is the first detailed 
and systematic numerical analysis on the influence of the 
oblateness coefficient on the character of orbits and this is 
exactly the novelty and the contribution of the current work. 

For several values of the oblateness coefficient in the last 
three Hill’s regions configurations we defined dense uniform 
grids of 1024 X 1024 initial conditions regularly distributed 
on the 0 < 0 part of the configuration ix,y) plane inside 
the area allowed by the value of the total orbital energy. 
All orbits were launched with initial conditions inside the 
scattering region, which in our case was a square grid with 
-2 < x,y < 2. For the numerical integration of the orbits in 
each type of grid, we needed about between 11 hours and 5 
days of CPU time on a Pentium Dual-Core 2.2 GHz PC, de¬ 
pending on the escape and collisional rates of orbits in each 
case. For each initial condition, the maximum time of the 
numerical integration was set to be equal to 10"^ time units 
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however, when a particle escaped or collided with one of 
the two primaries the numerical integration was effectively 
ended and proceeded to the next available initial condition. 

In this article we provide quantitative information re¬ 
garding the escape and collisional dynamics in the restricted 
three-body problem with oblateness. The main numerical re¬ 
sults of our research can be summarized as follows: 


1. It was observed that as the value of the oblateness coef¬ 
ficient increases the size of the stability islands around 
primary 1 is reduced and at relatively high values of the 
oblateness coefficient there is no indication of regular 
motion around the oblate primary, while at the same time 
the stability region around spherical primary 2 is also re¬ 
duced. 

2. The collisional basins to primary 1 were found to sig¬ 
nificantly increase in size with increasing value of the 
oblateness coefficient, while the basins formed by initial 
conditions of orbits which collide to the primary 2 are 
practically unaffected by the shifting on the value of the 
oblateness coefficient. 

3. It was detected that in the vicinity of the center of the 
oblate primary 1 a portion of orbits collide to the same 
primary almost immediately. Furthermore, the amount 
of these type of orbits increases with increasing oblate¬ 
ness. On the other hand this behaviour do not apply to 
the case of the spherical primary 2 . 

4. Our calculations reveal that as the primary body 1 be¬ 
comes more and more oblate the fractality of the planes 
is reduced and the boundaries between the escaping and 
collisional basins appear to become smoother. 

5. We presented evidence that the oblateness coefficient also 
influences the escape as well as the collision time of the 
orbits. In particular, both types of times are reduced as 
the value of the oblateness coefficient increases. 


Judging by the detailed and novel outcomes we may say 
that our task has been successfully completed. We hope that 
the present numerical analysis and the corresponding results 
to be useful in the field of escape dynamics in the restricted 
three-body problem with oblateness. The outcomes as well 
as the conclusions of the present research are considered, as 
an initial effort and also as a promising step in the task of un¬ 
derstanding the escape mechanism of orbits in this interest¬ 
ing version of the classical three-body problem. Taking into 
account that our results are encouraging, it is in our future 
plans to properly modify our dynamical model in order to 
expand our investigation into three dimensions and explore 
the entire six-dimensional phase thus revealing the influence 
of the oblateness coefficient on the orbital structure. 
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